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[i] We quantify gravity changes after great earthquakes present within the 10 year long 
time series of monthly Gravity Recovery and Climate Experiment (GRACE) gravity fields. 

Using spherical harmonic normal-mode formulation, the respective source parameters of 
moment tensor and double-couple were estimated. For the 2004 Sumatra-Andaman 
earthquake, the gravity data indicate a composite moment of 1 .2 x 1 0 23 N m with a dip of 
10°, in agreement with the estimate obtained at ultralong seismic periods. For the 2010 
Maule earthquake, the GRACE solutions range from 2.0 to 2.7 x 10 22 N m for dips of 
12°-24° and centroid depths within the lower crust. For the 2011 Tohoku-Oki earthquake, 
the estimated scalar moments range from 4.1 to 6.1 x 10 22 Nm, with dips of 9°-19° and 
centroid depths within the lower crust. For the 2012 Indian Ocean strike-slip earthquakes, 
the gravity data delineate a composite moment of 1.9 x 10 22 N m regardless of the centroid 
depth, comparing favorably with the total moment of the main ruptures and aftershocks. 

The smallest event we successfully analyzed with GRACE was the 2007 Bengkulu 
earthquake with M 0 ~ 5.0 x 10 21 N m. We found that the gravity data constrain the focal 
mechanism with the centroid only within the upper and lower crustal layers for thrust 
events. Deeper sources (i.e., in the upper mantle) could not reproduce the gravity 
observation as the larger rigidity and bulk modulus at mantle depths inhibit the interior 
from changing its volume, thus reducing the negative gravity component. Focal 
mechanisms and seismic moments obtained in this study represent the behavior of the 
sources on temporal and spatial scales exceeding the seismic and geodetic spectrum. 

Citation: Han, S.-C., R. Riva, J. Sauber, and E. Okal (2013), Source parameter inversion for recent great earthquakes from 
a decade-long observation of global gravity fields, J. Geophys. Res. Solid Earth , 118, 1240-1267, doi:10.1002/jgrb.501 16. 


1. Introduction 

[ 2 ] Large-scale processes of geophysical and climate- 
related mass redistribution cause changes in the gravitational 
potential field. By observing the relative motions of two 
identical satellites (orbiting proof masses), the Gravity 
Recovery and Climate Experiment (GRACE) mission has 
been mapping the spatial distribution of surface and interior 
mass flux and transport as well as adjustments in the Earth 
system since its launch in 2002 [Tap ley et al., 2005], As 
one of such processes, earthquakes cause variations in the 
gravitational potential field at a spatial scale up to some 
thousands of kilometers and at temporal scales of seconds 
to decades, by radiating seismic energy and deforming the 
surface and interior permanently and gradually. 
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[ 3 ] Traditional coseismic and postseismic observations 
have measured surface displacements such as ground, sea- 
floor, and sea-surface motions using a variety of instruments: 
seismometers, strainmeters, leveling, GPS, Interferometric 
Synthetic Aperture Radar (InSAR), seafloor transponders, 
and tsunami gauges. These have been the primary tools used 
to understand the rupture dynamics, the spatial extent of 
slip, and gradual postseismic changes. However, space- 
borne gravimetric observations are also sensitive, in partic- 
ular, to interior deformation in the broader region affected 
by the rupture, including the oceanic environment largely 
inaccessible by traditional measurements. Specifically, they 
could fill in the seldom-observed long wavelength spectrum 
of earthquake observations as a complement to surface 
geodetic measurements and seismic data. In addition, they 
reflect an average deformation over a time window much 
longer than that accessible from seismic data which is 
limited in principle by the period of the Earth’s gravest 
mode, 0 S 2 (3232 s). We examine this new type of earth- 
quake observations from GRACE gravimetry, which are 
sensitive to changes in gravitational potential which we 
express in the formalism of the Earth’s normal modes, as 
an alternative to the study of vertical and horizontal 
displacements, in order to seek additional information on 
earthquake mechanisms and a new perspective on earth- 
quake-related processes. 
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[ 4 ] We scrutinize the observations of global gravitational 
potential change after the following recent megathrust earth- 
quakes: 2004 Sumatra- Andaman (M w = 9.15 to 9.3), charac- 
terized by a slip greater than 3 m over more than 1200 km of 
fault length [e.g., Tsai etal., 2005; Chlieh etal., 2007]; 2007 
Bengkulu (Sumatra; M w = 8.5) [e.g., Borrero et al., 2009]; 
2010 Maule, Chile (M w = 8.8) with slip greater than 3 m over 
450 km [e.g., Pollitz, et al., 2011]; and 2011 Tohoku-Oki 
(M w « 9.0) with slip greater than 5 m over 330 km [e.g., 
Simons et al., 2011]. In addition to these thrust sources, and 
for the first time, we analyze gravitational perturbations due 
to the sequence of great strike-slip earthquakes (cumulative 
M w iv 8.7) which ruptured in 2012 100-200 km southwest of 
the Sumatra subduction zone [e.g., Yue et al., 2012], 

[ 5 ] Such gravity changes are caused dominantly by large- 
scale density change in the region surrounding the rupture 
and by surface deformation, including interaction with the 
ocean, as suggested by Han et al. [2006] and others there- 
after. The surface deformation, however, would be expected 
to be small in large-scale gravity change for vertical strike- 
slip earthquakes. Sea level changes due to coseismic gravita- 
tional potential changes have been addressed more recently 
by de Linage et al. [2009], Broerse et al. [2011], and 
Cambiotti et al. [2011]. While most of these analyses, in- 
cluding Heki and Matsuo [2010], Matsuo and Heki [2011], 
and Wang et al. [2012], used the monthly time series of 
GRACE gravity maps, Han et al. [2010, 2011] also directly 
exploited the fundamental observations of range-rate change 
after earthquakes to optimize spatial and temporal resolu- 
tion. Han et al. [2011] and Cambiotti et al. [2012] inverted 
the GRACE range-rate and monthly gravity data, respec- 
tively, to quantify the fault parameters of the 2011 
Tohoku-Oki earthquake. Various numerical modeling 
approaches from a half-space [Okubo, 1992] to a layered 
spherical Earth model [Pollitz, 1996], considering the ocean 
layer and sea level feedback [Broerse et al., 2011; Cambiotti 
et al., 2011], have been used to compute synthetic gravity 
changes. Furthermore, the postseismic gravity change after 
the 2004 Sumatra-Andaman rupture was also identified and 
such observations were interpreted considering the Earth’s 
rheological response to mega-earthquakes [Han et al., 
2008; Panet et al., 2010; Hoechner et al., 2011] and the 
diffusion of water in the mantle [Ogawa and Heki, 2007]. 

[6] In this study, we present a theoretical model of the 
coseismic gravitational potential for a spherical Earth on the 
basis of gravitational potential normal mode summation, as 
also envisaged by Chao and Gross [1987] and de Linage 
et al. [2009]. We particularly examine the gravitational char- 
acteristics of normal modes as a function of the Earth’s elastic 
structure and the earthquake source depth. We evaluate and 
discuss the effects of surficial layering (including the ocean) 
and density change on the gravity change. Next, we formulate 
inverse models to estimate the moment tensor components and 
subsequently fault parameters from the GRACE data for each 
earthquake. We use the readily available level 2 (L2) data 
products from the GRACE mission that implies monthly snap- 
shots of global gravitational potential changes with a spatial 
resolution of 400-500 km. Finally, we apply this inversion 
approach to characterize these earthquakes with (point) cen- 
troid source representation; discuss the fault solution estimates 
of moment, dip, strike, rake and depth, and the trade-offs 
among these parameters; and compare these great earthquakes 


based on homogeneous data, consistent gravity modeling, 
and uniform inversion methodology. Although the gradual 
postseismic changes after the rupture were also evident in 
the data, this study is focused on the coseismic change in the 
gravitational potential and leaves out the longer-term 
postseismic response to a future investigation. 

2. Coseismic Gravitational Potential Change Due 
to a Double-Couple Source 

[ 7 ] We derive the expression of coseismic changes in 
gravitational potential due to a point dislocation on the basis 
of gravitational potential nonnal mode summation. Our goal 
is to obtain the representation in terms of spherical harmonic 
coefficients for the coseismic gravitational potential changes 
as explicit functions of the parameters of a point-source 
double-couple (scalar seismic moment M 0 , dip 5, rake k, 
and strike <f>f ). In the end, we will use these results to analyze 
global gravity data such as from GRACE to invert earth- 
quake source processes. 

[s] The concept that gravity participates in the restoring 
force controlling the oscillation of a defonned elastic Earth 
was first expressed by Bromwich [1898], following Lamb's 
[1882] estimate of its period of free oscillation. Ever since 
Love's [1911] classical study, which provided the first theo- 
retical computation of the fundamental mode of a compress- 
ible, elastic gravitating Earth, all detailed computations of 
the Earth’s normal modes [e.g., Pekeris and Jarosch, 1958; 
Gilbert and Backus, 1968] have included the relevant varia- 
tions of its gravitational potential. Because they do not result 
in changes in its gravity field, the torsional modes of oscilla- 
tion of the Earth will be ignored from the rest of this paper. 

[ 9 ] The problem of the excitation of a nonnal mode by a 
seismic source was first described by Alterman et al. [1959] 
and given a simple and elegant fonnulation by Gilbert 
[1970]. This set the stage for the use of normal mode summa- 
tion to synthesize either seismic waves, which express 
the transient defonnation following an earthquake [e.g., 
Kanamori and Cipar, 1974] or geodetic displacements, which 
express permanent coseismic deformation [Pollitz, 1996], 
Thus, it should be possible to similarly describe any static or 
transient changes in the Earth’s gravitational potential on the 
basis of a summation of its normal modes. Indeed, this ap- 
proach goes back to Longman [1962, 1963] in the simpler case 
of the loading of the Earth by a point mass. 

[ 10 ] In the notation of Kanamori and Cipar [1974, 
Figure 10] and Kanamori and Given [1981, Figure 1], we 
recall that the radial displacement u r (as would be recorded 
by a vertical seismometer), excited by a point-source double- 
couple can be expanded as 

u r {r , 0, (j>; t) = {r)[-s R K 0 Pi <0 ( cos 6) + q R K x P h \ ( cos 9) 

n,l 

-PrK 2 Pi, 2(cos<9)] 

x [1 - cos{ n co,t)- exp(- n (0,t/2 n Q,)\. 

[ 11 ] The definition of the associated Legendre functions 
P !m used here is such that they are normalized to a constant 
integral of 4 ti on the sphere, which relates them to the 
Legendre functions Pf used, for example, by Kanamori 
and Cipar [1974] through 
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(a) m=0 (b) m=1 (c) m =2 



Figure 1. Spectrum of the three gravitational excitation functions evaluated at depths d s = 10, 20, and 
30km; (a) Fo(d s ) = ’S^Ko(d s )ys(a) (isotropic, m = 0), (b) F\(d s ) = E' (c4)y5 (a) (dipolar, m= 1), 

% n n 

and (c) F 2 (d s ) = } ^(dsjysja) (quadrupolar, m = 2). For Figure la, the isotropic function is also eval- 
uated at the upper and lower bounds of each distinct layer in PREM (i.e., upper crust, lower crust, and 
upper mantle) and shown as dashed lines. 


F '~ ~ \/* (2 ' +11 ! nw*’ <2) 

where k= 1 is for m = 0 and k=2 otherwise. The expansion 
in equation (1) uses a system of spherical harmonics whose 
pole is chosen at the seismic epicenter, and whose primary 
meridian is taken along the dislocation fault strike. Under 
this geometry, a point-source double-couple excites only 
the modes of azimuthal orders m = 0, ± 1, and ± 2, allowing 
the expansion ( 1 ) to contain only three terms inside the first 
bracket; s R , q R , and p R are then trigonometric parameters 
depending only on the geometry of the dislocation (dip angle 
5 and rake k), and on the longitude tfi of the receiver, mea- 
sured from the fault strike [Kanamori and Cipar, 1974, 
Figure 10]: 

sR = sin2 sin<5 cos<5, 

q R = sin2 cos2<5 sintf) + cos2 cosd coscj), 

Pr= cosz sind sin20 — sinl sin<5 cos<5 cos20, ' ' 

4> = 4>f- 4’r = <!>lon + <i>f-K, 

where (j)f is the azimuth of the fault strike, and <fi r the azi- 
muth of the small arc of great circle departing the source to 
the receiver, both measured clockwise from local north at 
the epicenter (alternatively, (p r = it — 0/ o „, where 0/ o „ would 
be the longitude of the receiver in a frame where the pole is 
at the epicenter, measured counterclockwise from the direc- 
tion of local south at the source). y\ is the vertical displace- 
ment component of the eigenfunction (dependent on degree l 
and overtone number /;) [Alterman et al., 1959; Saito, 1967]. 
The functions K 0 , K u and K 2 , of which explicit forms are 
given in Kanamori and Cipar [1974], describe the excitation 
of a full Earth oscillation by a seismic source, as such de- 
pends only on / and n as well as on the specific Earth model 
used, and of course on the depth d s of the source, but no 
longer on the particular orientation of either the source or 
the receiver. The temporal part expressed as the second 
bracket of equation (1) consists of a static term describing 
the permanent deformation of the Earth (expressed as the 
“1” in the bracket) and of a transient oscillatory term “seismic 
waves” corresponding to the attenuated sinusoid. An equiva- 
lent expression using moment tensor components can be 
found in Kanamori and Given [1981] and Pollitz [1996], 


[ 12 ] This formalism can be extended immediately to the 
case of the variations in gravitational potential ) 

by replacing the first component jq of the eigenfunction by 
the fifth one, y 5 , as introduced by Alterman et al. [1959], 
In the particular case of the permanent “static” variation of 
the potential at the Earth’s surface ( r — a ), one simply obtains 

0/ o „) = Mo-^2y5{a)-[-s R KoPi t o{ cos 6) (4) 

«,/ 

Fq R K\Pi \ ( cos 6) - PrK 2 P,, 2 { cos<9)]. 


[ 13 ] Note that y 5 is computed routinely as part of the solu- 
tion of the full elastogravitational eigenproblem of the free 
oscillation and that its inclusion is particularly critical to the 
accurate determination of the period of the gravest modes, in 
particular, the radial ones (/=0). However, the information 
contained in y 5 has rarely been used, since most applications 
of normal mode theory have been limited to seismology, 
before the advent of detailed large-scale gravity observations, 
notably with programs such as GRACE [Taplev et al., 2005], 

[ 14 ] By substituting equation (3), equation (4) can be re- 
written as follows: 


1 Ko, 4>i„n) = EX>™( cos(?){C/, m cos m(j) lon + 5/, m si nm<t> hn }, 

a 1 m = 0 

(5) 

where G is the universal gravitational constant, M the mass 
of the Earth, and a its radius. The five non-trivial dimension- 
less coefficients are defined, with the moment tensor compo- 
nents following the convention by Aki and Richards [1980] 
as follows: 


Ci. 0 — 


GM 


Cu = 


GM 


S,. 1 = 


Q, 2 = 


GM 


E Ko(d,)ys{a ) 

1 

J2Ki{d s )y 5 {a) 

n 

( d s)y5(a) 

n 

^ ^K 2 (d s )y 5 {a 


-M n 


GM 


( M , ) , 

f Mgg — M^ 

V 


(6a) 

(6b) 

(6c) 

(6d) 
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Si, 2 = 


a 

GM 


^ ~^K 2 {d s )y 5 {a ) 




where M rr =M 0 sin 2<5 sin X, 


(6e) 


M r g = —Mo ^ cos<5 cos2 cos 4>f + cos2<5 sin/, sinijyj , 

M r ^ = — M 0 cos<5 cos! sinijy + cos2<5 sinl cosily j , 

Mgg — M^ = —M 0 ^2 sint) cos/.- sin20y — sin2<5 sin2- cos2 </yJ , and 


Mgj, = —Mo ( - sin2<5 sin2 sin2<M + sin<5 cos2 cos2^ I . 


[is] In equation (6) and for each value of degree /, the sum- 
mations in the functions F m = K m (d s )ys («),(« = 0, 1,2) 

n 

are performed over the overtone number n on which the indi- 
vidual values ofy 5 (a) and the K m (d s ) will of course depend. In 
particular, these functions of / are now characteristic of the 
(laterally homogeneous) Earth model, but no longer depend 
on any property of the seismic source (excluding depth 
through the K m (d s )), and as such, they can be precomputed 
for any combination of / and d s . 

[16] These functions are closely related to the coefficients 

lm(S ) 

y x used to expand radial deformation by Pollitz [1996, 
equation (3a)], except that we express them with sines and 
cosines as opposed to imaginary exponentials and that 
Pollitz [1996] considered only elastic (and not gravitational) 
restoring forces. Pollitz [1992] showed in considerable detail 
(including the effect of gravity) that it is possible to compute 
directly the functions F m from the equations of equilibrium 
of the deformed Earth, rather than explicitly summing over 
an a priori infinite number of overtones n, as in equation 
(6). We follow here an approach similar to Piersanti 
et al.’s [1995] and Pollitz et al.’s [2006] studies, which con- 
sists of casting the static problem of the equilibrium of a lay- 
ered gravitating elastic Earth in the same fonnalism as the 
dynamic one for a free oscillation, but of course in the ab- 
sence of any time-dependent terms, following in the foot- 
steps of Longman [1963] and Smylie and Mansinha [1971]. 

[ 17 ] In this framework and given a point-source double- 
couple of prescribed geometry and amplitude ( M 0 ), imbedded 
at the depth d s in an elastic Earth, Pollitz ’ [1996] formal- 
ism consists of expanding its field of discontinuities 
(obtained under a representation theorem) onto spherical 
harmonics, and of defining, for each degree l and orders 
m = 0, 1, and 2, an equivalent discontinuity at the depth d s in 
the four components of the eigenvector y of the static problem 
(which remains independent of m). In the presence of gravity, 
we generalize it by imposing the continuity of the additional 
components, y 5 and for r=r s ( r s =a — d s ), thus imposing 
six boundary conditions on the full six-dimensional 
elastogravitational eigenvector y. Following Pollitz et al. 
[2006] and because y is continuous at all depths d ^ d s , except 
at the core-mantle boundary as discussed for example by 
Smylie and Mansinha [1971], it can be integrated continu- 
ously upwards to r = r~ from the core-mantle boundary, 
where initial conditions feature three degrees of freedom, 
detailed for example by Vermeersen et al. [1996], in the 
framework of Chinnerv’s [1975] explanation of the so-called 
Longman’s paradox. Similarly, it can be integrated from the 
surface of the Earth (where there are again three degrees of 
freedom) downwards to r+. For each m = 0, 1, and 2, the six 


boundary conditions at r=r s define the proper combination 
of initial conditions from which the value of 75 ( 0 ) can be com- 
puted and equated to F m . 

[is] As mentioned above, equations in (6) are relatively 
simple because they are written in a particular system of 
spherical hannonics, where the pole is at the seismic epicen- 
ter and the prime meridian oriented along the local south, 
since under this geometry, a point-source double-couple ex- 
cites only azimuthal orders m = 0, ± 1, and ± 2. That frame 
will hereafter be identified with a superscript s. However, 
when studying the gravity field of the Earth, the convention 
is to use the geographic spherical harmonics, defined about 
its axis of rotation (the North Pole) and the prime Greenwich 
meridian, hereafter be identified with a superscript g. There- 
fore, one has to project one system of spherical harmonics 
onto the other, in order to express equation ( 6 ) in the geo- 
graphic harmonics; this relatively classical problem [e.g., 
Sato, 1950; Brink and Satchler, 1968; Stein and Geller, 
1977] is carried out as a succession of two solid rotations: 
The first one (of amplitude 0 S , the colatitude of the epicenter) 
is taken about the pole of the geographic meridian going 
through the epicenter and brings back the pole of the spher- 
ical harmonics to the geographic North Pole, and the second 
one (of amplitude — (j) s , the opposite of the longitude of the 
epicenter) is taken about the axis of rotation of the Earth 
and restores the Greenwich meridian as the primary one. 

[ 19 ] Regrouping the coefficients C; and S \ , in equation (6) in 
the fonn of a vector Xj (Xf _ 2 = S ,, = S,, u X[ Q = C/, 0 , 
Xf x = C/,i, and Xf 2 = C/, 2 ), the projection onto the geographic 
harmonics is expressed as a new vector Xf : 

2 

xf m (m 0 , <j) f 5, X, d s , e s , 0,) = J2 E U e °’ W; (M), 4ft, X, ds) 

-l<m<l, (7) 

where the factor E expresses the combination of two rotations. 
Note that despite its apparent complexity, equation (7) remains 
a completely linear operation on the various components of 
the vector Xj. 

[ 20 ] A fundamental aspect of equation (7) is that while the 
potential field resulting from a dislocation expressed as a 
point-source double-couple could be expanded on hannonics 
of orders m = 0, ± 1, and ± 2 in the source-based system 
(X) being five-dimensional), it will project on all orders 
(— l<m<l) in the geographic system centered at the North 
Pole (Xf being (21+ l)-dimensional). Note that a similar ap- 
proach is used when studying the splitting of the free oscil- 
lations of the Earth due to rotation and ellipticity, with the 
identical result that normal modes excited by an earthquake 
are split into all of their (2 /+ 1) singlets (— / < m < l) [Stein 
and Geller, 1977]. 

3. Coseismic Gravitational Response Functions 

[ 21 ] We characterize the coseismic changes in gravita- 
tional potential 1 j/ through the functions 

F m(d s ) = ^ ^K m (d s )y 5 (a ), ( 8 ) 

n 

with orders m = 0, 1 , 2 for isotropic, dipolar, and quadrupolar 
excitations, respectively. They represent the Earth’s gravi- 
metric response to faulting by a double-couple and depend 
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only on Earth structure, degree /, and source depth d s . Figure 1 
presents examples of these functions, computed for the Pre- 
liminary Reference Earth Model (PREM) model [Dziewonski 
and Anderson, 1981] at three representative depths (10, 20, 
and 30 km), sampling the upper crustal (3—15 km), lower 
crustal (15-24 km), and upper mantle (24-80 km) layers of 
the model. 

[ 22 ] As shown in Figure 1, the behavior with / of the three 
functions is significantly different. Both F\ and F 2 are 
weakly dependent on depth and monotonic (respectively, in- 
creasing and decreasing functions of /) functions, while F 0 
strongly depends on depth, to the extent that it even reverses 
sign, with the nodal degree decreasing with increasing 
depth. This strong dependence of F 0 on depth is further in- 
vestigated by computing its values at the upper and lower 
bounds of each source layer, shown as dashed lines on 
Figure la; they suggest a moderate dependence of F 0 with 
depth within a layer, but a significant one between layers, 
which illustrates the critical dependence of the coefficients 
K 0 on the elastic moduli at the source. The spatial pattern 
of F m at different depths is shown in Figure 2, truncated at 
degree 50 (i.e., spatial resolution of 400 km). The isotropic 
term, F 0 , is dominant when the source is at the lower depth, 
while the non-isotropics, F, and F 2 , become more pro- 
nounced as the source depth increases. 

[ 23 ] Several features of Figures 1 and 2 are straightfor- 
ward to interpret in the context of classical results of normal 
mode seismology. In particular, we recall that the strike-slip 
coefficient K 2 is proportional to y 2 /r s . In the vicinity of the 
Earth’s surface, where the shear stress _y 4 vanishes, the 


logarithmic derivative of this quantity takes the value 

n ■ 1? (t) = ~ ~ - it’ on the order of “ l/a for 1 > 

and for branches of modes featuring a reasonably circular 
ground motion at the surface (e.g., fundamental Rayleigh 
waves for which this ratio approaches 1.5). Thus, the charac- 
teristic depth for the decay of K 2 near the surface is on the 
order of all, which means that K 2 is essentially independent 
of depth for shallow earthquakes, a well-known property of 
mantle Rayleigh waves. This explains the lack of depen- 
dence of F 2 on depth, as exhibited in Figure lc. By contrast, 
the coefficient K\, proportional to J/ 4 // 2 , vanishes at the sur- 
face, leading to the classical singularity in excitation by 
dip-slip components at shallow depths [e.g., Kanamori and 
Given, 1981]. For shallow sources ( d s much less than one 
wavelength), y 4 is expected to be continuous and grow line- 
arly with depth; however, will suffer discontinuities in the 
various layers of the Earth, leading to a much slower growth 
(or even an irregular behavior) of K x with depth. As for the 
isotropic coefficient K 0 , involved (together with K 2 ) in the 
case of a dipping (e.g., thrust) fault, its more complex ex- 
pression as a function of the eigenvector y precludes any 
simple interpretation of its variation with depth. 

[ 24 ] The physical origin of the change in gravitational po- 
tential at the surface, _v s (a), can be separated into two contri- 
butions, namely (i) a displacement of the free boundary of 
the Earth, resulting in the filling of an initially void volume 
by material of finite density and (ii) the effect on i/i at the sur- 
face of changes in density within the interior of the planet, 
resulting from its deformation, including that of surfaces of 
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Figure 2. Spatial patterns of the functions F„,(d s ) in the vicinity of the source, evaluated at depths of 10, 
20, and 30 km. 
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discontinuity, such as the crust-mantle interface. For (ii), the 
gravitational effect of the (coseismic) change of discontinu- 
ity in interior density stratification is considerably smaller 
than the interior density change by compression and/or dila- 
tation [Pollitz, 1997; Han et al., 2006], In the limit of large- 
scale gravity changes (a few 100 km), the effect from (ii) is 
as large as (i), while (i) overwhelms (ii) at much smaller 
scales, as shown in Han et al. [2006], 

[ 25 ] The change in gravitational potential due to the defor- 
mation of the surface, F , is in the nature of a Bouguer 
effect, and can be easily computed from the radial dis- 
placement component of the eigenvector, yfa). In the limit 
yx < all (i.e., the vertical deformation being smaller than its 
lateral dimension), the contributing material can be modeled 
as a thin layer of thickness y\ and density A p (crustal density 
for inland earthquakes, reduced by seawater density for 
undersea ones), leading the classical result [Jeffreys, 1976; 
Turcotte et al., 1981]: 

F rn ) = (2/+^) - yi W ’ m = °’ 1 ’ 2 - (9) 

[ 26 ] The degree-dependent scale factor multiplying the y\ 
eigenfunction is needed to convert the thin layer mass anom- 
aly A py\(a) into the gravitational potential anomaly at the 
surface. The contribution from the interior deformation 
(mostly due to coseismic dilatation or compression of the 
material surrounding the source) is computed by removing 
this superficial (or Bouguer) contribution from the full vari- 
ation F m as 

pi 1 ) — p _ p( B ) 

= Y^K m (d s ) (y s (a) - (a)j , m = 0, 1,2. (10) 

[ 27 ] Although we did include self-gravitation and loading 
when computing the functions F m and for the sake of simplic- 
ity, they are not discussed here, as their contributions are con- 
siderably smaller than the simple Bouguer effect expressed by 
equation (10); see de Linage et al. [2009] and Cambiotti et al. 
[2011] for details. The degree-dependent behaviors of 
equations (9) and (10) at three different depths are shown in 
Figure 3, where the sign for the surface deformation function 
for m = 0 was reversed in the plot. As for the case of total grav- 
ity change in Figure 1, the depth-dependence is found to be 


significant in the isotropic case, where the interior deformation 
shows even greater depth-dependence. For the isotropic com- 
ponent, which is directly related to dilatation at the source (see 

Appendix A), the interior defonnation Fq \ (m = 0) contrib- 
utes more signal than the Bouguer term F ^ when the source 
locates at shallower depth (10 km, within the upper crust) 
where the material is more compressible. In contrast, Fq B) be- 
comes prominent when the source is deeper (i.e., 30 km, 
within the upper mantle) where the material is less compress- 
ible. In the dipolar case (m=l), the change in potential is 
largely due to the Bouguer term F {B} at all depths, while in 
the quadrupolar case (m = 2), F^ remains prominent at all 
depths. 

[ 28 ] We computed the coseismic potential change from a 
synthetic double-couple at different depths (using a com- 
pressible Earth model with an ocean layer). A centroid with 
a scalar moment M 0 = 5.0 x 10 22 Nm (M w = 9.0) and strike 
cj)/=180°, dip §=15°, and rake /, = 90' was used. The 
gravity changes (up to degree and order 40, or equivalently 
a spatial resolution of 500 km) from the same seismic source 
but at different depths are compared in Figure 4. The effect 
of the isotropic term ( m = 0), mostly responsible for a large 
central negative anomaly, decreases with increasing depths 
and, as a result, the other terms (m= 1,2) that are relatively 
constant with depth become more prominent in the total 
coseismic gravity change. 

3.1. Moment Dip Trade-off 

[ 29 ] It has been long known in the seismological commu- 
nity that the inversion of seismic moment tensors suffers a 
singularity for shallow sources [e.g., Kanamori and Given, 
1981], This is due to the fact that the dipolar coefficient 
K], proportional to the shear stress j> 4 , vanishes at the surface 
with the result that the two components M rip and M rS do not 
excite seismic modes at or near the surface. Conversely, they 
cannot be resolved for a surficial source (or are poorly re- 
solved for a shallow source) from seismograms. This classi- 
cal singularity results in a trade-off between scalar moment 
M 0 and dip angle <5. 

[ 30 ] For thrust earthquakes (Aw 90°), the observations of 
isotropic (m = 0) and quadrupolar (m = 2) gravity response 
(equations (6a) and (6d)) constrain only M 0 sin2<5. The 
complementary information required to separate the scalar 


(a) m=0 (b) m=1 (c) m =2 



Degree 

Figure 3. Same as Figure 1, with function F m (d s ) separated into their surface (F'F ) and interior (F'F ) 
contributions to the total gravitational potential. Note that the sign of Fq B) (shown as solid line) has been 
reversed in Figure 3 a for plotting clarity. 
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Figure 4. Examples of synthetic gravity changes from a double-couple seismic source at different depths 
with the layered Earth model are shown in the top panels. The focal mechanism parameters are given by 
M 0 = 5.0 x 10 22 Nm, (j)f= 180°, <5=15°, and 1 = 90°. The source is located at the center of each diagram 
within the upper crust (10 km), lower crust (20 km), and upper mantle (30 km). The lower panels show 
the contribution to the total gravity change Figures 4a-4c from each order (m = 0, 1, and 2, respectively). 


moment M 0 and dip 5 comes from the observations of the 
dipolar (m= 1 ) gravity response (equation ( 6 c)) containing 
a factor of M 0 cos2<5. The resolvability of M 0 and 8 from 
(noisy) gravity measurements would, therefore, be depen- 
dent on the ratio of S u (and C /4 depending on strike) 
to Ci ft. Taking a reasonable range for shallow-dip events 
(5° <<5<30°), the ratio is approximately given by 


O 



( 


o 


\ 


^ATi(^)y5(a) 


^-KoK).V 5 (>) 


\ 


, which is expected to 


go to zero with all K\ coefficients at the surface. Assuming 
that the accuracy of satellite gravity measurements depends 


only on degree /, the signal-to-noise ratio (SNR) of 5/4 and 
C/ o can be related through: 


SNR{5/4 } _ 

SNR{C/ i0 } ~ ~ F 0 (d s ,l) ' 


[ 31 ] Figure 5 presents the ratios (11) for various source 
depths. As expected and by the analogy with seismological 
inversions [Kanamori and Given , 1981], for very shallow 
sources, this ratio remains small over the entire bandwidth 
(/ < 50) accessible by the GRACE satellites. At such shallow 
depths the dipolar term (m= 1) is simply too poorly excited 
to resolve both M 0 and <5 from a noisy large-scale satellite 
gravity dataset. Flowever, for other sources within the lower 
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Figure 5. Ratio of the dipolar gravity response F\ (d s ) to 
the isotropic response F 0 (d s ) as function of degree 1. This ra- 
tio characterizes the resolvability of moment from dip angle 
when inverting gravity observations, as a function of degree 
l for various source depths. 

crust and the upper mantle, SNR {.S'/j } becomes comparable 
to SNR{C /j0 }, likely allowing the resolution of M 0 and <5 
from the GRACE gravity observations. 

4. Observations of Global Gravity Field and 
Spatial Localization 

[ 32 ] As part of the GRACE project, one month’s worth of 
orbital tracking data of the GRACE satellites have been rou- 
tinely analyzed to obtain an “average” snapshot of the global 
gravity field at each month since April 2002. Such a snapshot 
is represented by and provided as a set of spherical harmonic 
coefficients, known as L2 product. In this study, we will exam- 
ine those L2 spherical harmonic coefficient data from April 
2002 to September 2012. The L2 data have been continuously 
updated with the improved processing strategy and back- 
ground models including atmosphere and ocean variability. 
We used the most recent Release-5 (RL5) L2 data product 
except for 2002, for which the Release-4 (RL4) product was 
used, since the RL5 product has, at present, not yet been 
processed for these years [Bettadpur et al., 2012], 

[ 33 ] The coseismic gravitational perturbation tends to be 
spatially localized around the epicenter, and thus, its energy 
is smeared out among all spherical harmonic coefficients. 
Due to other gravitational variations from continental-scale 
mass variations (mainly in association with seasonal climate 
changes), the earthquake signal may not be apparent directly 
from the L2 data. Such "local” signal can be better delin- 
eated from the “global” spherical harmonic coefficients by 
applying a spatial window around the region of interest. 
Simons et al. [1997] and Wieczorek and Simons [2005] 
developed an optimal windowing function to extract the spa- 
tially confined signal by minimizing the spatial truncation 
(leakage) effect. Flan and Simons [2008] applied this tech- 
nique successfully to identify the 2004 Sumatra-Andaman 
earthquake signal out of the spherical harmonic coefficients 
from GRACE. Han and Ditmar [2008] discussed how 
the SNR of non-stationary signals may be substantially 
underestimated, if not localized, and suggested a better 


way to quantify the SNR of the local signal from the spher- 
ical harmonic coefficient data under stationary noise. 

[ 34 ] The optimal spatial windowing function h(6), often iso- 
tropic and band-limited, is found by maximizing its energy 
within a confined region of interest [Wieczorek and Simons, 
2005]. Once the windowing function is determined, one can 
obtain spherical harmonic coefficients of the spatially local- 
ized signals using a “convolution-like” formula provided in 
equation (10) of Wieczorek and Simons [2005] and given as: 

Cl = (-If E E hj C im V(2/+l)(2/+l)(2/+l) 

J=° HH\ (12) 

i j l\ ( i j l \ 

0 0 0 ) \m 0 —m J ’ 

where L h is the maximum expansion degree of the 
windowing function and hj is the expansion coefficients of 
the zonal (isotropic) window h(9). The last two parentheses 
are Wigner-3j functions. C im is the (original) spherical 
harmonic coefficient and C\ m is the localized coefficient 
highlighting the signals over the region. 

[ 35 ] As implied by equation (12), the coefficient of the 
localized (or windowed) signal is nothing but a linear combi- 
nation of the original coefficients with the weights defined 
by the choice of the windowing function. The localized coef- 
ficient of degree / is computed with the original coefficients 
within the bandwidth [max(0, / — L h ), (l+L h )]. If the original 
signal is known only to a certain degree such as L s , it can be 
readily seen from the upper bound of the summation in 
equation (12) that the permissible range of the localized coef- 
ficients is limited to L s — L h . Also the windowed spectra of the 
low degrees (/ <L h ) can be biased as noted in Section 5.1 of 
Wieczorek and Simons [2005]. Therefore, we use the localized 
spherical harmonic coefficients within the bandwidth, L h + 1 
to L s — L h , from monthly GRACE L2 data. 

[36] Figure 6 shows signals for recent great earthquakes 
(the 2004 Sumatra-Andaman, 2007 Bengkulu, 2010 Maule 
(Chile), 2011 Tohoku-Oki, and 2012 Indian Ocean events) 
that are in a detectable range with GRACE gravity data. 
The spatial maps of coseismic gravity change computed 
from Global Centroid Moment Tensor (GCMT) solutions 
of each respective event are shown in Figures 6a-6e, with 
the spherical circle that delineates the area of localization 
with the spherical cap of radius 0 h . The cross section of each 
spatial windowing function is shown in Figures 6f-6j, 
respectively. The relative amplitude is presented over the 
spherical angular distance 9 from the center of the cap. 
The expansion degree L h is fixed at 20 and the cap radius 
9 h is chosen to capture most of the (expected) coseismic 
gravity signal, mainly depending on the size of the earthquake. 
This provides the only optimally localized windowing func- 
tion for each earthquake (that is, the Shannon number is equal 
to 2 as in Wieczorek and Simons [2005]). The amplitude of 
each windowing function becomes approximately 1% of the 
maximum amplitude beyond the spherical cap of radius 9 h , 
indicating that h(9) is nearly perfectly concentrated within 
the spherical cap. The concentration ratio y, showing how 
much power of the window function is concentrated within 
the spherical cap, approaches 99% for all cases. 

[ 37 ] The temporal variability (in terms of root-mean-square, 
RMS) of the monthly GRACE L2 data, after localization is 
applied in the respective areas, is shown in Figures 6k-6o, 
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Figure 6. (a-e) Synthetic gravity changes computed from centroid moment tensor (CMT) solutions for the 
2004 Sumatra-Andaman, 2007 Bengkulu, 2010 Maule, 2011 Tohoku-Oki, and 2012 Indian Ocean earth- 
quakes, respectively. The black circle delineates the spherical cap of radius 9 \ as the region of localization 
where the spatial windowing function li{9) is concentrated; (f-j) cross sections showing the relative amplitude 
of li(9), expanded to degree 20 ( L h = 20), but with a different cap radius ( 9 h ) for each event; (k-o) degree RMS 
spectrum (square root of power spectrum in pGal) of synthetic coseismic gravitational potential (blue line) 
and of average variability of monthly GRACE L2 data (red line) after applying an identical spatial windowing 
function in each region. The data variability (red line) includes the gravity signal variability (e.g., seasonal 
change) as well as the GRACE instrument noise. The ratio of the (predicted) coseismic gravitational pertur- 
bation to the GRACE data variability approaches 6, 0.5, 1, 3, and 1.5, respectively for each event. 


for each earthquake, within the bandwidth, L h + 1 to L s — L h 
(i.e., from 21 to 40 in our case since L s = 60 and L h = 20). 
This is suggestive of inherent temporal variability of the 
gravity field (mass redistribution) in each region and of the 
observational noise. The gravity signal strength predicted 
from GCMT solution of each event is also depicted after 
the same localization is applied. The ratio of the (predicted) 
coseismic gravitational perturbation to the GRACE data 
variability approaches 6, 0.5, 1, 3, and 1.5, respectively 
for each event, within most of the bandwidth from degrees 
21 to 40. In Section 6, we will examine the time series of 
these localized L2 coefficient observations and invert them 
to determine the fault parameters of each earthquake. 

5. Linear and Nonlinear Inversion 

[38] So far, we examined the forward model of the gravi- 
tational potential change in response to a point-source 


double-couple, expressed in terms of the moment tensor as 
in equation (6) and in particular, its dependence on depth. 
The spatial localization of global GRACE gravitational 
potential data was introduced for identifying and analyzing 
coseismic gravity changes around various earthquake 
regions. In this section, we discuss how we determine the 
moment tensor from the “localized” GRACE data and, sub- 
sequently, invert the double-couple source parameters from 
the moment tensor solutions. 

5.1. Linear Estimation of Moment Tensor from Gravity 
Data 

[ 39 ] The estimation of the moment tensor components 
from a series of geopotential observations is straightforward 
due to its linearity, when the data are expressed in terms of a 
series of spherical hannonic coefficients, as discussed in 
Section 2. In order to invert the dataset of localized GRACE 
gravitational potential coefficients into the moment tensor m. 
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we express the forward problem as a succession of three 
operations. 

[ 40 ] First, as expressed in equation (6), we multiply the 
tensor of Green’s functions evaluated at the seismic source 
depth d s by the moment tensor: 

X s = F(c/ V )m, (13) 

where F(c/ V ) regroups, for all used values of /, the five 
functions F m discussed in Section 3; m = ( M rr , M,. g , M rq b , M gg 
— M g 0 ) represents the five independent unknown mo- 
ment tensor components; the resulting vector X' s regroups, 
for all values of /, the spherical harmonic coefficients in 
equation (6). Second, we apply the rotations expressing the 
change of frames, and obtain 

X'aRfiWX 1 , (14) 

where the dimension of the vector X" is the sum of (21 + 1) 
over the range [Imin, Imax] used in the inversion, and R 
(0 s ,4>s) is a rotation matrix consisting of all the elements 
E l mi of equation (7). Third, the spatial localization discussed 
in Section 4 is applied, leading to 

X £ »L(0„0,)X*, (15) 

where L,(O s ,(f> s ) is a matrix expressing the operation in 

equation (12). 

[ 41 ] At this point, X L regroups the gravitational potential 
coefficients expressed in the geographic frame, and localized 
around the epicenter (9 s ,(t> s ), as excited by the double-couple 
source described by m. The inversion then consists of using 
the GRACE dataset of observed coseismic change in gravi- 
tational potential x (itself expended onto spherical har- 
monics) to solve for m through 

L(6> s , </> s )x = L(0„ 4> S )K(6 S , , 4> s )T(d s )m + e, (16) 

where e is the noise in the GRACE dataset. Note that L ap- 
pears on both sides of (16) because e is determined only 
when the localized data Lx is sought from noisy time series 
of GRACE data. Finally, simple least squares is used to de- 
termine the estimate of the moment tensor components, m, 
and its covariance matrix, C{m}. 

5.2. Nonlinear Estimation of Double-Couple From 
Moment Tensor 

[ 42 ] The determination of double-couple source parameters 
(i.e., M 0 , (j>f, 5, and X) from moment tensor components con- 
sists of solving backwards the equations in (6). While point- 
source double-couples are characterized by four parameters, 
they do not constitute a four-dimensional vector space, and 
moment tensor inversions have to be carried out in a five- 
dimensional vector space, fostering both nonlinearity and 
non-uniqueness. We use a standard approach in seismology 
to derive two nodal plane solutions of strike, dip, and rake 
by using a computer code (bb.m, a Matlab code written by 
Oliver Boyd, based on mij2d.f, a FORTRAN code by Chen Ji) 
available from the website [www.ceri.memphis.edu/people/ 
oldboyd/Softwarc/Softwarc.html |. They are used as initial so- 
lutions for our refined inversion considering variable uncer- 
tainties in moment tensor components as described below. 

[ 43 ] First of all, we examine correlation among the fault 
parameters. For low-dip (8 « 0) earthquakes as in most 


of cases in this study, the moment tensor components are 
approximated (to the first order) to M„. « M 0 2S sin 2, 
M r0 « - M 0 cos (<j>f- X), M r<t> « M 0 sm(cj) f — X), M gg - « 

—M(,28 sin(20y— X), and M g< j , ra — Mq8 cos(2 (f>f— X). It indi- 
cates M,. g and M r<t) would be dominant and (j)f— X (not individ- 
ually) would be better-constrained from the moment tensor. 
The parameters such as (f)f and X would be tightly coupled in 
a manner that (j)j — X is constant. Han et al. [201 1] discussed 
this trade-off between strike and rake found from various 
seismic solutions and GRACE data inversion for the 2011 
Tohoku-Oki earthquake. 

[ 44 ] The correlation among double-couple parameters can 
be quantified by examining a covariance matrix of the 
double-couple solution. We derive it by introducing a small 
perturbation in the double-couple parameters as follows: 


M 0 = M 0 (l + drf), (17a) 

<t>/ = 4>f + d<j>f, (17b) 

8 = S + dS, (17c) 

2 = 2 + dX. (17d) 


[ 45 ] The perturbation in the vector of moment tensor com- 
ponents Am due to perturbation in double-couple parameters 
d is simply written as Am « Bd, where B is a matrix 

consisting of partial derivatives, e.g., evalu- 

ated for nominal values of Mo, <5, X, and <pf. The Am includes, 
e.g., AM,.,., andA(Mflg — M^), and finally d consists of dr], 
dd, dX, and difif. Therefore, the covariance matrix of the 

double-couple parameters, C{d}, is computed from the co- 
variance matrix of the moment tensor solution, C{m}, by 
error propagation: 


C 




(18) 


[46] The correlation matrix is simply computed by re- 

A 

scaling all elements of C{d} by p t • = c,y / y/cijcj] , where 
Cjj is a component of the z-th row and j-th column of the 
matrix. This is a metric we examine using the initial fault 
plane solutions from moment tensor estimates for all earth- 
quakes in this study. 

[ 47 ] The (iterative) least square refinement to the initial 


solution of Mo, <5, X, and <j)f, (available from the eigenvalue 
and eigenvector decomposition of the moment tensor matrix, 

A A „ 

C{d}B r 


i.e., from the code bb.m), is found by d 
C{m } 1_1 


(ih — m) where m is a moment tensor component 


vector computed with the initial double-couple parameters. 
This solution takes care of the variable uncertainties and 
correlation in the moment tensor estimates from GRACE. 
Depending on the correlation structure among four param- 
eters computed with equation (18), this solution may need 
to be constrained. We will elaborate on this procedure for 
each of the earthquakes in Section 6. 

[48] We do not attempt to solve for the centroid location 
(9 s ,(j) s ) and depth d s because their dependence on gravity 
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data, as expressed in equation (16), is algebraically complex. 
However, we test various depths d s , and at each depth, we 
solve for the corresponding moment tensor and the best 
double couple. We adopt the centroid location ( 0 S ,4 > S ) deter- 
mined from other seismic CMT or finite fault solutions. 

6. Inversion for Moment Tensor and Fault 
Parameters 

6.1. 2004 Sumatra-Andaman Earthquake 

[ 49 ] We analyzed the entire time series of monthly global 
gravity fields from April 2002 to September 2012 (in tenns 
of spherical harmonic coefficients) in order to separate, in 
the data obtained after the rupture, the gravity signal due 
to the earthquake from the background temporal variation in 
the gravity potential. We fit the time series using the mean, 
annual and semiannual sinusoids, and a Heaviside step 
(for coseismic) and a logarithmic function (for postseismic) 

simultaneously. The logarithmic term, log f 1 + 150 days ) 

where t is the time elapsed since the rupture, is used to elim- 
inate the effect of viscous mantle relaxation on a time scale 
significantly longer than that expected from a seismic 
source, even allowing for slow coseismic components. The 
linear and quadratic components were not included in the re- 
gression because they are correlated with the postseismic 
trends; however, the coseismic step estimates are not 
affected, regardless of whether the linear and quadratic com- 
ponents are included or not (by virtue of sufficiently long 
time series). In the tune series regression, we did not use 
the data in 2012 to avoid any potential influence by the 
nearby 2012 strike-slip earthquakes. 

[ 50 ] Figure 7 shows the time series of the GRACE L2 co- 
efficients, localized within the domain shown in Figure 6a. 
The localized coefficient at degree / is computed with the 
neighboring coefficients, and thus creates correlation among 
the localized coefficients over degree (but independent over 
order), effectively reducing random noise in the original L2 
data. For example, the localized coefficient at degree 30 was 
influenced by the original GRACE coefficients at degrees 10 
to 50, because we used the windowing function expanded 
with the maximum degree 20 ( L h = 20). The (localized) coef- 
ficients were plotted in the (epicentral) coordinate system 
where the z axis locates at the center of the spatial window 
function (5.4°N, 93.8°E). It is the location where the largest 
moment was released according to the finite fault model 
used in Han et al. [2006], In this rotated coordinate system, 
the gravity coefficients at orders 0, 1, and 2 are directly 
related to the moment tensor components of the centroid at 
the pole (i.e., 5.4°N, 93.8°E in a geographic coordinate 
system), as described in equation (6). Almost all coefficients 
between degrees 21 and 40 present large perturbations due 
to the earthquake in episodic change as well as gradual 
change afterwards. Although postseismic observations last 
long (>7 years), they need to be carefully examined since 
they might be affected by background gravity changes due 
to inter-annual climate variability. The longer postseismic 
time series reveal that they may be better modeled with a 
logarithmic function than the exponential function used in 
Han et al. [2008] for the analysis of only 2.5 years of 
postseismic data. The coseismic SNR is found highest from 
the coefficients of C/j, Si\, and S^ 2 , for M r0 , M r ,f„ and M g 


respectively, and degrades for C/ >2 and C t (h related to the diag- 
onal components of the moment tensor, M rr and M gg — 
respectively. 

[ 51 ] The coseismic step and its error estimates for each of 
the localized coefficients were used for (linear) inversion of 
the moment tensor components m and its covariance matrix, 
cjmj. From the condition M rr + M gg + = 0, we found 

three diagonal components. The GRACE estimates at depth 
of 20 km and GCMT are compared in Figure 8, indicating 
that GCMT provides smaller estimates of the components, 
particularly M r( /,. The corresponding focal mechanism dia- 
gram (i.e., “beachball”) was drawn by finding two fault 
planes using the program mentioned in Section 5.1. We 
use them only as initial solutions for our inversion of the 
double-couple parameters reflecting the error characteristics 
of the moment tensor estimates from GRACE. 

[ 52 ] First of all, the correlation matrix as discussed in 
Section 5.2 was computed for both fault-plane solutions and 
shown in Figure 9. As suspected from the low dip fault plane, 
rake and strike are tightly coupled with nearly unity correla- 
tion, while the other fault plane solution does not yield such 
strong coupling between them. This pair of correlation charac- 
teristics is typical for other thrust events considered in this 
study (note, however, that the correlation matrix looks differ- 
ent for the strike-slip event that will be discussed later). 

[ 53 ] Due to strong correlation between strike and rake, we 
fixed strike <f>f a priori and solved for M 0 , 8, and k simulta- 
neously from the moment tensor and its covariance estimates 
following the procedure in Section 5.2. The convergence 
was always obtained within not more than three iterations. 
The solutions at various strikes for both fault planes are 
shown in Figure 10. By changing (/yin 1° steps around the 
initial strike, we found that double-couple solutions for the 
primary fault plane with 320° <(/y<360° fit equally well 
the GRACE observations. The rake parameter changes line- 
arly with strike, while the other parameters of M 0 and <5 are 
relatively constant. The trade-off between strike and rake is 
only found for the plane with low dip angle. For the conju- 
gate fault plane, the secondary double-couple solution can 
be delineated with </y= 150°±5°, M o ~750 x 10 2O Nm, 
(5-82°, and 1-88°. 

[ 54 ] We tested sensitivity of the fault solutions to depth. In 
this case, we fixed the strike at 340° and changed depths 
within the lower and upper crust and upper mantle (i.e., from 
10 to 30 km). The solutions in the crustal layers yielded var- 
iance reduction (VR) = 0.9 1-0.95, while for solutions in the 
upper mantle (not shown), VR= 0.70-0.75 was substantially 
lower (Figure 11). Any solution within the crustal layer fit 
GRACE observations equally well with the ones in the 
lower crustal layer being even better. With increasing 
depths, 8 gradually increases, while M 0 decreases and k is 
practically unchanged. The solutions at shallower depths, 
particularly within the upper crust, show larger changes in 
M 0 with the variation in 8 or depth than do the deeper 
sources, as expressed by a steeper slope in the d s — M 0 or 
8 — M 0 plot in Figure 1 1 . This indicates that M 0 is not robustly 
resolved for the upper crustal solutions. As expected, for 
shallow thrust earthquakes, M 0 sin2(5 is better-constrained 
[Kanamori and Given, 1981; Tsai et al., 2011]. As depth 
increases (i.e., lower crust), however, the effect becomes less 
acute, and the resolution of M 0 and 8 improves. 
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Figure 7. Monthly time series of the GRACE L2 data after applying the spatial localization over the 
Sumatra-Andaman earthquake region (solid blue line) and the seasonal and inter-seasonal fit (solid red 
line). Arbitrary offsets were added for clarity. The data residual (black line with error bars) was computed 
by subtracting the fit (red line) from the data (blue line). The associated error bar was estimated from the fit 
obtained between 2002 and 2011. The data residuals were subsequently analyzed using the Heaviside step 
and logarithmic functions for delineation of coseismic and postseismic changes, respectively (solid green 
line). The GRACE coefficients were rotated to the epicentral coordinate system where the z axis locates at 
5.4°N and 93.8°E and the five lowest order (m = 0, 1, and 2; for both C/„, and S/,„) coefficients are shown at 
different degrees between 21 and 40. From the top to bottom, these coefficients related to the moment 
tensor components, M m M r0 , M,.^, M 00 — and M g respectively. 
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Mrt (M r0 ), Mrp (M,.,;,), and Mtp (Mg^) from GRACE and GCMT. (right) For GRACE, the centroid depth 
was 20 km which is the focal mechanism corresponding to two possible fault planes (NP1 and NP2). This 
initial fault solution was computed from GRACE estimates of moment tensor without consideration of 
variable uncertainties in different tensor components. 
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Figure 9. The correlation matrices of the fault (double-couple) parameter estimates for the two fault planes. 
Note nearly unity correlation between strike and rake for the fault plane (NPl) with a low dip angle. 


[55] The single CMT solutions from GRACE at various 
depths within the crustal layers show scalar moments vary- 
ing from 6.4 to 10 x 10 22 Nm and dip angles varying from 
3° to 13°. This is a considerably higher estimate than the 
GCMT (4 x 10 22 Nm and a dip of 8°). In terms of amplitude 
and spatial pattern of gravity change (Figure 12), the GRACE 
solutions at various depths are all consistent, with differences 
smaller than 2 pGal at the spatial resolution of 500 km. The 
primary negative anomaly appears at the back-arc region and 
the secondary positive anomaly appears offshore. 

[56] As a trial of multiple-source solutions, we introduced 
four centroids for the GRACE data analysis. We first deter- 
mined the centroid locations corresponding to the major asper- 
ities based on a finite fault model [Han et al., 2006]. The fault 
parameters (M 0 , <5, and /.) of each of the four centroids at a 
depth of 20 km were simultaneously estimated with strike 
fixed as 340° for all centroids and the dip angle was assumed 
to be uniform for the all centroids (that is, total nine indepen- 
dent parameters were estimated: four scalar moments, four 


rakes, and one dip for four centroids). Expectedly, the new 
solution with four sources shows improvement, yielding 
VR=0.99. For the case of the multiple CMT solutions from 
GRACE, the total moment is 1 160 x 10 2O Nm and the dip is 
1 0° , producing a gravity change of similar (or slightly larger) 
amplitude as the single centroid solution, but elongated along 
strike with the rakes decreasing northward (Figure 13). The 
overall distribution of the moments with changing geometry 
(combination of strike and rake) in N-S direction is consistent 
with the five moment tensor solutions found from long-period 
seismic data [Tsai et al., 2005]. 

6.2. 2007 Bengkulu Earthquake 

[ 57 ] This earthquake of moment 5 x 10 21 Nm (M w = 8.4) 
[Borrero et al., 2009] is the smallest one that was detected 
by GRACE gravity data during the period from 2002 to 
2012. The GCMT reports this event as a shallow-dip thrust 
rupture, striking (j)f= 327° at depth of 23.3 km. It released 
the largest moment (~40 x 10 2O Nm) in the component of 
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Figure 10. The fault plane solutions determined from the moment tensor and its error covariance esti- 
mates for both planes with various strikes. The depth was fixed as 20 km for all cases. In each case, the 
strike was fixed a priori and Mo, <5, and a were estimated simultaneously. The solutions with strike ranging 
between 320° and 360° are acceptable with VR>0.94 for the primary fault plane and the ones with strike 
from 145° to 155° are acceptable for the secondary fault plane. Note that 2 is linearly dependent of the 
strike while M 0 and (5 are relatively consistent for the primary fault plane, as also predicted from the cor- 
relation matrix (Figure 9). 



Figure 11. The fault plane solution at various depths within the crustal layers. The strike was fixed at 
340° for all cases. In each case, depth was fixed a priori and the moment tensor was determined at the 
corresponding depth and then the double-couple parameters were estimated from it under the constraint 
4>f= 340°. The solutions with the centroid below the crustal layer do not fit the GRACE data. Mo and 8 
are highly dependent on depth while 2 is not. 


M rS . The 10 years of GRACE gravity time series localized 
over the region (Figure 6b) identified a statistically signifi- 
cant change in C/,i coefficients, before and after the earth- 
quake epoch as shown in Figure 14. Such coseismic jump 
in the coefficients estimates the M r/) moment release of 
4.8 ± 0.3 x 10 21 N m at a depth of 24 km, favorably compar- 
ing with the GCMT moment estimate. The second largest 
moment release (~2.5 x 10 21 N m) is expected from the 


M,.0 component; however, noisier S^\ coefficients of the 
GRACE gravity data make it difficult to identify this subtle 
change. GRACE satellites measure the gravity change in 
the N-S direction (controlled by M r( j) much better than in 
the E-W direction (controlled by M, y; ,) due to their better 
sampling of along track orbit perturbation [ Han et al., 2005]. 

[58] The 2005 Nias earthquake ruptured 3 months after the 
2004 Sumatra-Andaman earthquake and had a moment 
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Figure 12. (a) Surface gravity changes computed from the 
GRACE solution for the 2004 Sumatra-Andaman earth- 
quake, at a spatial resolution of 500 km (or spherical har- 
monic degree up to 40). The particular solution with the 
centroid depth of 20 km was used, (b-d) Difference in sur- 
face gravity between the solution at 20 km depth and the 
ones at 10, 15, and 24 km, respectively. The estimates of 
scalar moment, strike, dip, and rake for each solution are 
presented at the bottom of each panel. The location of the 
centroid is 5.4°N and 93.8°E. 

magnitude twice as large as the 2007 event (1 10 x 10 2O N m). 
However, this event was difficult to separate in the monthly 
GRACE data after the large coseismic and postseismic gravity 
change of the 2004 Sumatra-Andaman earthquake. 

6.3. 2010 Maule (Chile) Earthquake 

[59] A localization of the GRACE L2 time series in the 
Chilean region (Figure 6c) reveals coseismic and (possibly) 
postseismic gravity changes in association with the 2010 
Maule earthquake (Figure 15). It effectively diminishes 
other signals and highlights the gravity variations confined 
to the specific region, including seasonal and inter-annual 
trends. This seasonal and inter-annual gravity change was 
removed as in the previous case. The residuals were then an- 
alyzed with the coseismic step and postseismic logarithmic 
functions. Likewise, the estimates of the coseismic step 
and its error from the localized coefficients within the degree 
band (/ = 2 1 — 40) were used for the moment tensor and fault 
parameter inversion. 

[60] The estimated moment tensor components from 
GRACE are compared with other seismic solutions in 
Figure 16. For this earthquake, the moment release in M r $ 
(or S /,1 coefficients) is much larger than M r0 (or C/_ 1 coeffi- 
cients), indicating an overall strike along N-S. Compared 
to seismic solutions, GRACE overestimates M n . and M gg . 
The fault plane solutions and focal mechanisms are depicted 
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Figure 13. Same as Figure 12a for the 2004 Sumatra- 
Andaman earthquake multiple centroid solutions at a depth 
of 20 km. Total gravity change from four centroid moment 
tensors estimated from GRACE. The strike was fixed at 
340° for all centroids and the dip is assumed to be the same 
for all centroids. The individual location of the centroids 
(shown as red star) were determined by approximating the 
finite fault model used in Han et al. [2006], They are (3.4°N, 
94.9°E), (5.4°N, 93.8°E), (8.4°N, 92.7°E), and (11.8°N, 
92.5°E). Note that the overall distribution of the scalar 
moment and rake (or strike) northward is consistent with the 
long-period seismic solutions found in Tsai et al. [2005], 

on Figure 16, using the moment tensor estimates at a depth 
of 20 km. As in the previous case, these initial double-couple 
(fault) solutions were derived without considering variable 
uncertainties in the moment estimates. We also found that 
the correlation matrix for the low-dip fault plane (primary 
one) shows a tight coupling between strike and rake while 
the conjugate plane with larger dip does not, just like for 
the 2004 Sumatra-Andaman earthquake. 

[61] Our final double-couple solutions considering the error 
covariance of the moment tensor estimates were determined 
by nonlinear inversion starting with the initial solutions. 
Figure 1 7 shows various solutions for two fault planes at depth 
of 20 km with variable strikes that are fixed a priori prior 
to each inversion. If we delineate the solutions yielding 
VR>0.88, M 0 varies between 215 and 222xl0 2O Nm 
(deceasing with strike), <5 varies between 17° and 20° (increas- 
ing with strike), and rake changes linearly between 80° and 
115° (increasing with strike), when (/y changes from 5° to 
35°. The conjugate fault solution can be described with 
M o = 220±7 x 10 20 Nm, <5 = 72°±2°, 2 = 85°±3°, when 
(j)/= 185° ± 7°. 

[62] The GRACE solutions at depths in the upper mantle 
(25-30 km) indicate poor agreement with the GRACE 
observations, yielding only VR < 40%; on the other hand, 
the inversions in the upper and lower crust (10-24 km) yield 
VR~90%. The GRACE solutions obtained within the 
crustal layers are presented in Figure 18. As in the previous 
case, Mq and <5 change significantly with depth while X is 
more or less constant. The depths of the Maule earthquake 
seismic solutions range from 24 to 35 km. The crustal 
thickness from CRUST 2.0 [Bassin et al., 2000] over the 
Maule area is 32 km, thicker than the PREM model we 
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Figure 14. Same as Figure 7 for the 2007 Bengkulu earthquake. Note that the C / , coefficients present a 
significant step before and after the epoch of the event in 2007. This particular earthquake was strongest in 
the component of M rB . 


used. Therefore, the GRACE solutions at the lower bound 
of the lower crust would be preferable. This particular 


solution is also consistent with most of other solutions in 
terms of the dip angle and moment [Lay et al., 2010; 
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Figure 15. Same as Figure 7 for the 2010 Maule (Chile) earthquake. The GRACE coefficients were ro- 
tated in the coordinate system where the z axis locates at 35.6°S and 286. 9°E. 


G Shao et al., Preliminary Slip Model of the Feb 27, 2010 
Mw 8.9 Maule, Chile Earthquake, 2010, http://www.geol.ucsb. 
edu/faculty/ji/big_earthquakes/20 1 0/02/27/chile_2_27.html; 


G. P. Hayes 2010, Finite Fault Model, Updated Result of 
the Feb 27, 2010 Mw 8.8 Maule, Chile Earthquake, 
http://earthquake.usgs.gov/earthquakes/eqinthenews/20 1 0/ 
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Figure 16. Same as Figure 8 for the 2010 Maule (Chile) earthquake. 



Figure 17. Same as Figure 10 for the 2010 Maule (Chile) earthquake. The depth was fixed as 20 km for 
all cases. The solutions with strike ranging from 5 to 35° are acceptable, with VR>0.88 for the primary 
fault plane and the ones with strike from 182° to 192° are acceptable for the secondary fault plane. 


us2010tfan/finite_fault.php]. Also, the GRACE solutions 
show M 0 consistently larger than for other long-period 
point-source CMT solutions, agreeing better with the finite 
fault models. The GRACE solutions within the upper crust 
fit the observations equally well; however, the estimate of M 0 


sin 2d is smaller than that of the solutions within the lower 
crust presumably due to the influence of a lower rigidity. 
The GRACE solutions within the upper crust layer, although 
presented with large variations in M 0 , are more robust in terms 
ofM 0 sin2d. 
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Figure 18. Same as Figure 11 for the 20 10 Maule (Chile) earthquake. The strike was fixed at 19° for all cases. 


[63] The gravity change was computed using the GRACE 
CMT solution at a depth of 24 km (Figure 19). It shows, 
likewise, the primary, negative, anomaly on land (back-arc 
region) and secondary, positive, anomaly offshore. The 
differences in terms of gravity (at 500 km resolution) among 
the solutions at depths of 10, 15, and 20 km are less than 
1 pGal, while the dip angle estimates vary from 6 to 24° 
and the moment estimates from 207 to 276 x 10 20 Nm, 
depending on depth. 

[64] For this event, we note that the lateral surface density 
heterogeneity due to the South American continent (being 
different from the uniform ocean layer used in our calcula- 
tion) might be important. We made some ad hoc computa- 
tions to roughly quantify the effect of land by combining 
the ocean anomaly computed with the uniform ocean layer 
and the land anomaly computed without the ocean layer. 
The eigenfunctions were computed using two different ID 
Earth models with and without the ocean layer. We found 
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Figure 19. Same as Figure 12 for the 2010 Maule (Chile) 
earthquake. The location of the centroid is 35.6°S and 
286. 9°E. 


the lateral heterogeneity may cause a difference in the 
negative gravity on land by 2 pGal at maximum. We defer 
more reasonable assessment of the surface heterogeneity 
in gravity to the computation of eigenfunctions with a 3D 
Earth model and sea level equation such as in Broerse 
et al. [2011], 

6.4. 2011 Tohoku-Oki (Japan) Earthquake 

[65] Figure 20 presents the time series of the same 
GRACE L2 coefficients, but localized around the area 
shown in Figure 6d, at m = 0,1, and 2 and the regression fits. 
After removing the mean, linear, and quadratic polynomials 
and annual and semi-annual sinusoids, the residual data are 
analyzed for the coseismic and postseismic changes. All five 
moment tensor components present significant coseismic 
and postseismic signature in the time series. Compared to 
the time series around the Maule earthquake, this area is 
affected less by systematic variations in gravity. 

[66] The coseismic step and its error estimates for each of 
the localized coefficients were used for inversion of the 
moment tensor components at various depths. The solutions 
at 20 km are shown in Figure 21 along with other seismic 
solutions. The focal mechanism parameters were obtained 
for the primary and conjugate fault planes as also shown in 
Figure 21. All moment tensor components from GRACE at 
20 km are in good agreement with seismic solutions. Once 
again, the fault plane solution with low dip angle yields nearly 
unity correlation between rake and strike, while it is not 
the case for the other fault plane. Therefore, we fix the 
strike a priori to obtain the double-couple solutions, and 
then test various strikes (Figure 22). Solutions featuring 
Af 0 = 450^-70 x 10 2O N m (increasing with strike), 5 = 14° -5° 
(decreasing with strike), X changing linearly from 65° to 90° 
(increasing with strike), and <py= 175°-95° all fit the GRACE 
observations equally well, with VR> 98%. The conjugate fault 
solution can be described with M 0 = 450±20 x 10 2O Nm, 
<5 = 77° ± 1°, X = 95° ± 3°, when </y= 1 8° ± 3°. 

[67] The solutions with the centroid in the upper mantle 
(depth greater than 24 km in PREM) indicated consistently 
poorer agreement with the GRACE observations, yielding 
VR ~ 50%, compared to the ones within the upper and lower 
crust with VR > 90%. Therefore, we ruled out the solutions 
within the upper mantle and only show in Figure 23 the 
other solutions located within the crustal layers. The upper 
crustal solutions show dip angles of 7° and less and scalar 
moments of 3 80-600 x 10 2CI N m. The solutions within the 
lower crust indicate a dip ranging from 9° to 19° and M 0 
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Figure 20. Same as Figure 7 for the 2011 Tohoku-Oki (Japan) earthquake. The GRACE coefficients 
were rotated in the coordinate system where the z axis locates at 38.5°N and 142.6°E. Efigher noise in 
2002 is due to the use of RL04 GRACE products. The quality of the product is improved with RL05 as 
shown in the rest of the years. 


ranging 410 to 620 x 10 2o Nm, depending on the precise 
depth. In both cases, the shallower the depth in each layer, 
the smaller the dip and the larger the moment. When com- 
pared with other seismic and GPS-based solutions including 
long-period CMT and the moment-weighted CMT from var- 
ious finite fault models [Shao et al., 2011; Simons et al., 
2011; Hayes, 2011; Ammon et al., 2011], only the GRACE 
solutions within the lower crust are consistent with those 


alternative solutions in moment and dip, as well as depth 
(except for the USGS CMT in depth). 

[68] The coseismic gravity change (at a spatial resolution 
of 500 km), computed from the GRACE solution for a depth 
of 20 km, is shown in Figure 24a. It locates the primary neg- 
ative anomaly in the back-arc region and the secondary pos- 
itive anomaly (with one third the amplitude of the primary 
anomaly) near the trench. The alternative GRACE solutions 
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Figure 21. Same as Figure 8 for the 201 1 Tohoku-Oki (Japan) earthquake. 



Figure 22. Same as Figure 10 for the 2011 Tohoku-Oki (Japan) earthquake. The depth was fixed at 
20 km for all cases. The solutions with strike ranging from 175° to 195° are acceptable with VR>0.98 for 
the primary fault plane and the ones with strike from 15° to 23° are acceptable for the secondary fault plane. 
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Figure 23. Same as Figure 1 1 for the 2011 Tohoku-Oki (Japan) earthquake. The strike was fixed at 190° 
for all cases. 
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Figure 24. Same as Figure 12 for the 2011 Tohoku-Oki 
(Japan) earthquake. The location of the centroid is 38.5°N 
and 142. 6°E. 


at 10 (middle of the upper crust), 15 (upper bound of the 
lower crust), and 24 km (lower bound of the lower crust) 
are represented on Figures 24b-24d as deviations from the 
field of Figure 24a. They all remain within 1-2 pGal of the 
solution at 20 km. 


6.5. 2012 Indian Ocean Strike-Slip Earthquake 

[ 69 ] We examined the same GRACE coefficients for the 
same period but localized around the area where a sequence 
of great earthquakes ruptured in the Indian Ocean off 
Northern Sumatra on 11 April 2012 (Figure 6e). Figure 25 
shows the time series from 2006 to September 2012. It 
was found that the long-term trend particularly evident in 
the dipolar coefficients, C u and .S’/ 1 , is primarily due to the 
on-going viscoelastic deformation after the 2004 Sumatra- 
Andaman earthquake (e.g., numerical modeling such as in 
Han et al. [2008]). Unlike the previous earthquakes, the di- 
polar coefficients (for M r0 and M,.^) do not show any signif- 
icant coseismic change, while the quadratic coefficients such 
as M g< j, and {M ge — M^) present a significant coseismic sig- 
nature, albeit the latter component being noisier. 

[ 70 ] The GRACE moment tensor estimates at depth 30 km 
were presented with other seismic solutions of the M w = 8.6 
main rupture on 1 1 April in Figure 26. It is clear that the 
GRACE estimates, indicating a composite of the ruptures, 
are substantially larger than seismic CMT solutions. The 
M g $ estimate indicates a moment twice as large as derived 
from conventional seismic solutions. The initial focal mech- 
anism is also computed and shown in Figure 26. The corre- 
lation matrix for the double-couple parameters does not 
show any strong coupling, for either fault plane, as shown 


in Figure 27. This is due to the large dip angle of both fault 
planes. 

[ 71 ] Therefore, we inverted four double-couple parameters 
simultaneously from the moment tensor and its covariance 
estimates. Figure 28 shows these solutions for the primary 
and secondary fault planes at various depths from 10 to 
50 km. All of them fit the GRACE data equally well, with 
VR>0.95. For the primary fault plane, M 0 = 165-210 x 10 2O Nm 
(decreasing with depth), (5 = 78°-88° (increasing with depth), 
X = 1 83°-190° (decreasing with depth), and </>y is more or less 
constant at 1 1 1° — 1 13°, as depth changes from 10 to 50km. 
For the secondary fault plane, within the same depth range, 
Mq and 8 are similar as for the first fault plane, while X varies 
from —2° to —14° (decreasing with depth) and <j>f is more or 
less constant at 19°-22°. 

[ 72 ] These double-couple solutions are relatively consis- 
tent regardless of depth, which is distinctly different from 
the previous thrust events. It can be explained as follows: 
The gravity change due to a strike-slip earthquake is charac- 
terized mostly by the quadratic gravity coefficients such as 
Q 2 and 5/, 2 . Those coefficients are excited through the func- 
tion K 2 in Kanamori and Cipar [1974], As we extensively 
examined in Section 3, this function is nearly independent 
of depth at the range we considered in this study. These 
GRACE solutions are consistent with, or slightly larger 
than, the cumulative moment release during the first 2 h after 
the main rupture (139 x 10 20 Nm [Yue, et al., 2012]) com- 
bined with that of the aftershock (40 x 10 20 Nm [Duputel 
et al., 2012]). 

[ 73 ] The coseismic gravity change (at a spatial resolution 
of 500 km) computed from the GRACE solution at a depth 
of 30 km (WNW-ESE plane) is shown in Figure 29a. It 
locates the primary quadrupolar anomaly at the epicenter. 
Although the free-air anomaly (surface and interior deforma- 
tion) is shown in the figure, the Bouguer gravity anomaly 
(interior deformation) is very similar since the surface verti- 
cal deformation is small, as can be expected from Figure 3c. 
It indicates that the gravity anomaly from this earthquake is 
mostly from the interior process of compression (yielding 
positive gravity) and dilatation (yielding negative gravity). 
The alternative GRACE solutions at 30 km (a conjugate 
NNE-SSW plane), 10 km (WNW-ESE plane), and 50 km 
(WNW-ESE plane) are shown on Figures 29b-29d as devi- 
ations from the field of Figure 29a, respectively. The differ- 
ences among them remain less than 1 pGal, meaning that all 
these double-couple solutions explain the GRACE gravity 
change equally well. 

7. Summary and Discussion 

[ 74 ] By analogy with the representation of seismic waves 
using normal mode summation, we obtained an analytic form 
for the coseismic gravitational potential changes induced by a 
point-source double-couple on the basis of the gravitational 
potential component of the elastogravitational normal mode 
eigenfunctions. A double-couple source excites isotropic, 
dipolar, and quadrupolar patterns of gravitational change that 
are expressed in terms of spherical harmonics (centered at 
the epicenter) in orders 0, 1, and 2, respectively. Each spheri- 
cal harmonic coefficient is related to the moment tensor 
components of M rr , M rS , M,.^, M gg — M^, and M^, in a sim- 
ple linear manner. 
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Figure 25. Same as Figure 7 for the 2012 Indian Ocean earthquake. The GRACE coefficients were 
rotated in the coordinate system where the z axis locates at 3.8°N and 101. 0°E. The time series was drawn 
to start from 2006 to highlight 4 months of the GRACE observations after the rupture. 


[ 75 ] Theoretical computations of the gravitational poten- 
tial Green’s functions F m (d s ) for shallow sources within lay- 
ered Earth models yield the following results: 

[ 76 ] 1. The fault depth (related to rigidity and bulk modu- 
lus) is critical to characterizing the spatial pattern and ampli- 
tude of the coseismic gravity perturbation, mostly through 
the isotropic component. For example, isotropic gravity re- 
sponses from sources within the crustal layers and the upper 


mantle layer can be opposite in sign. In contrast, the 
quadrupolar (“strike-slip”) function is independent of depth. 

[ 77 ] 2. The interior deformation (density change) is more 
sensitive to the source depth and is larger than the surface 
deformation (Bouguer effect) for the shallower sources 
(d s < 20 km). 

[78] 3. For shallow thrust sources located in the upper 
crust, large-scale gravity data constrain only M 0 sm28. 
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Figure 26. Same as Figure 8 for the 2012 Indian Ocean earthquake. Note that GRACE indicates a 
cumulative moment release from all events while various CMT solutions involve only the major event. 
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Figure 27. Same as Figure 9 for the 2012 Indian Ocean earthquake. Unlike the previous thrust events 
characterized by a shallow-dipping fault plane, there is no tight coupling among the double-couple param- 
eters for both fault planes. 


[ 79 ] 4. For low-dip thrust events, strike and rake are 
tightly coupled in a manner such that only their difference 
is constrained. 

[so] All of these indicate that there might be trade-offs 
among the source parameters: scalar moment and dip with 
fault depth and strike and rake when dip is small. 

[si] We developed an inversion method to estimate the 
fault parameters from global gravitational potential observa- 
tions expressed as spherical harmonic coefficients, and 
applied it to the 2004 Sumatra- Andaman, 2010 Maule (Chile), 
2011 Tohoku-Oki (Japan), and 2012 Indian Ocean earth- 
quakes. A significant change in certain gravity coefficients 
related to M rS was also observed after the 2007 Bengkulu 
earthquake. For the thrust events, we found that centroid 
solutions within the upper and lower catst (d s < 25 km) give 
considerably better fits to the GRACE data than the deeper 
ones do located within the upper mantle (d s > 25 km). We 
interpret this result as due to an increase in bulk modulus 
(incompressibility) with depth responsible for the discordance 
in the resulting synthetics, as also found by Cambiotti et al. 
[2011], In addition, crustal solutions (both upper and lower) 
yield more robust values of the parameter P~ 1 Mq sin2<5, where 
\i s is the source rigidity, and thus of A u sin 28 (where A u is the 
seismic slip). Note that the latter is directly representative of 


the vertical component of slip in the limit angles 8 for which 
cos 8 is stationary. Other seismic and geodetic solutions are 
consistent with the GRACE solutions within the lower crust. 
For the strike-slip events in 2012, the double-couple solutions 
present no significant depth-dependency and no trade-off 
among the fault parameters. 

[ 82 ] In terms of spatial patterns of coseismic gravity 
changes, the GRACE observations consistently found pri- 
mary negative gravity anomalies in the back-arc region and 
smaller positive anomalies offshore, for the three thrust 
events in 2004, 2010, and 2011. The negative anomaly orig- 
inates mostly from the isotropic term expressing the influ- 
ence of dilatation through the excitation function K {] of 
Kanamori and Cipar [1974]. The deeper sources, located 
within the upper mantle with larger bulk and rigidity moduli 
than in the crust, cannot reproduce the larger negative grav- 
ity changes. In this respect, the satellite gravity observations 
from the megathrust events can be interpreted as expressing 
large-scale interior deformation associated with density 
changes (dilatation) inside the crustal layers. From the 
strike-slip event in 2012, the gravity change is, again, mostly 
from the interior process since the gravity effect from the 
vertical deformation at the crust-ocean interface is an order 
of magnitude smaller. The quadrupolar gravity pattern, 
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Figure 28. Same as Figure 11 for the 2012 Indian Ocean earthquake. In this earthquake, the simulta- 
neous solutions of the double-couple parameters were obtained at various depths for both fault planes 
(due to no tight coupling among the parameters). All the solutions produce VR>0.95. A distinct differ- 
ence from the previous thrust events is that the double-couple solutions of this strike-slip earthquake 
are relatively consistent regardless of depths. The scalar moment estimate from GRACE is consistent with 
the cumulative moment release from a series of the main ruptures and aftershocks. 


being positive in ENE and WSW quadrants and negative in 
NNW and SSE quadrants, indicates compression and dilata- 
tion, respectively, radially inside the Earth. 

[83] Our work also offers the possibility of investigating 
the presence of any ultraslow components in the spectrum 
of the seismic source. The largest ever earthquake recorded, 
the 1960 Chilean event was documented by Kanamori and 
Cipar [1974], Kanamori and Anderson [1975] and Cifuentes 
and Silver [1989] to include a slow precursor, which started 
15 min before the high-frequency rupture, and may have had a 
comparable moment. Similarly, a number of studies of the 
2004 Sumatra earthquake (including the modeling of its nor- 
mal modes) have shown that its very-long period moment, 
1.15 x 10 23 Nm, was more than 2.5 times greater than that 
derived from the standard GCMT algorithm, indicating an 
element of slowness also expressed in other characteristics of 
its source [Stein and Okal, 2005; Tsai et al., 2005; Choy and 
Boatwright, 2007]. Preliminary evidence also suggests that 
the 1 964 Alaska earthquake, and possibly the Rat Island earth- 
quake of 1965, may have featured a similar component of 
source slowness [Nettles et al., 2005]. By contrast, the normal 
modes of the 2010 Chilean and 2011 Tohoku events do not 
reveal any such slowness [Okal et al., 2012; Okal, 2012]. 

[84] Because the inversion of GRACE L2 data uses 
successive orbits for a month, the seismic moments resolved 
in the present study represent averages over a time window 
much longer than accessible from seismic data. The latter 


are limited in principle by the period of the Earth’s gravest 
mode, qS? (3232 s), and in practice most classical CMT in- 
versions are carried at periods not exceeding 300 s. Conse- 
quently, the comparison between our results and classical 
seismological ones can shed some light on the behavior of 
the source on time scales transgressing the seismic spectrum: 
a significant disparity between a GRACE moment and a 
(seismic) CMT solution would be a proxy for the existence 
of a slow source component releasing moment at frequencies 
below that of 0 S 2 , while obviously it could not resolve the 
details of such a component (e.g., whether it was precursory, 
coseismic or took place 10 days after the main shock). 

[85] 1. For the 2004 Sumatra-Andaman earthquake, the 
composite GRACE solution consisting of four sources elon- 
gated along strike yields a composite moment of 1 .2 x 1 0 23 
N m with a dip of 1 0° , in excellent agreement with the value 
of 1.2 x 10 23 Nm obtained at ultralong seismic periods by 
Tsai et al. [2005] and Okal and Stein [2009], and about 
2.5 to 3 times greater than the classical GCMT solution. 
We conclude that GRACE confirms the slow components 
occasionally documented at frequencies lower than those 
of CMT inversions, but fails to unearth even slower ones, 
with characteristic times of hours to days, whose existence 
could not have been a priori excluded. 

[86] 2. For the 2010 Maule, Chile earthquake, our solu- 
tions range from 2.0 to 2.7 x 10 22 Nm for dips between 
12° and 24°, and centroid depths between 15 and 24 km. 
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Figure 29. Same as Figure 12 for the 2012 Indian Ocean 
earthquake. The location of the centroid is 3.8°N and 
101.0°E. Figure 29 a stands for the WNW-ESE fault plane 
while Figure 29b the NNE-SSW plane at the same depth of 
30 km. Note that the surface gravity is not significantly differ- 
ent even with a wide range of variation in centroid depth. 


slips with a finite dimension of fault would be more suitable 
for such large earthquakes. Our quantification of the centroid 
sources from gravity observations compares, in general, fa- 
vorably with other seismic centroid solutions. When jointly 
analyzed with other terrestrial data (like GPS and seismic re- 
cords), GRACE gravity data will provide seldom-observed, 
broad-scale constraints on the spatial distribution of slip. 

[ 91 ] The gradual postseismic changes over a time frame 
of several years observed in most cases should not be 
overlooked. Continuous observations of gravitational poten- 
tial changes will be useful to document the rheological 
response of the Earth to great earthquakes and to advance 
our understanding of post-earthquake stress and strain redis- 
tribution. These data could provide, for example, an insight 
into bulk behavior: time-dependent (inelastic) bulk modulus 
[Rundle, 1978; Rundle and Smith, 1982] versus no signifi- 
cant bulk relaxation [Cohen, 1980a, 1980b; Ivins and 
Sammis, 1996; Pollitz, 1997]. These could be indispensable 
for understanding the large-scale behavior of the Earth’s 
shallow material, and eventually for earthquake hazard 
assessment and mitigation applications. Finally, the future 
GRACE follow-on mission equipped with enhanced instru- 
mentation [Watkins, 2011] may allow us to exploit gravita- 
tional potential data in the analysis of smaller, and thus more 
frequent seismic events than the ones studied here. The tran- 
sient gravitational perturbations available from LIB orbit 
perturbation data could deliver additional understanding at 
the seismic frequencies. 

Appendix A: On the Behavior of Isotropic Gravity 
Response From Shallow Sources 


These numbers compare favorably with finite fault models 
such as Lay et alls [2010] (M 0 = 2.1 x 10 22 Nm), and are 
only about 20% greater than the GCMT solution. 

[ 87 ] 3. In the case of the 2011 Tohoku-Oki earthquake, 
scalar moments estimated from GRACE data range from 
4.1 to 6.1 x 10 22 Nm, with dips between 9° and 19° and 
depths between 1 5 and 24 km. This is in excellent agreement 
with the published GCMT value (5.3 x 10 22 N m) and sup- 
ports OkaV s [2012] conclusion that the earthquake does 
not feature anomalous source slowness. 

[88] 4. For the 2012 Indian Ocean strike-slip earth- 
quakes, GRACE estimates the scalar moment at 1.9 x 
10 22 Nm at a depth of 30 km, 60% larger than the com- 
bined GCMT values for the mainshock and aftershock 
(0.9 x 10 22 Nm and 0.3 x 10 22 Nm respectively), but com- 
paring favorably with the total moment of 1.4 x 10 22 Nm 
of the finite fault solution of the main shock sequence by 
Yue et al. [2012] combined with the estimate of 0.3 x 10 22 N 
m for the aftershock [Duputel et al., 2012], 

[ 89 ] 5. The detection of the 2007 Bengkulu earthquake 
(M f) = 50 x 10 2O Nm; M w = 8.5) may demonstrate the lower 
bound for the use of GRACE gravity data to resolve earth- 
quake deformation. This was feasible due to the highest sen- 
sitivity of GRACE data to gravity change in N-S direction, 
where the spatial resolution could improve 400-500 km to 
as little as 200-300 km. 

[ 90 ] In this work, we attempted to characterize five great 
earthquakes in terns of a centroid (point) representation, 
although the laterally and vertically distributed moments or 


[ 92 ] In order to investigate the isotropic gravity response 
function F 0 (d s ), we examine the excitation function n K 0 , 
given in Kanamori and Cipar [1974] using the eigenfunction 
components defined by Alterman et al. [1959]. This can be re- 
written by introducing the dilatation component as (we specify 
here the degree / and overtone n, explicitly): 


n K oM) = 


21+1 


47i„(rj[/i + 1(1 + \)h\ 


{^,(40 -3,^(4)}, (Al) 


where ,,07 is the eigenfrequency; n X t is the dilatation 
eigenfunction that can be written as „X, = ’’’ 1 2 


; fi+i) 1 

'~T~ n y 3 .r ancl 

-2 X 


l nYl,l % r nYl.l 

_y w is the radial derivative of n y u : n y u = 
1 ,, 1 1 ) „ with the Lame 


(i + %Y r " yu + ( x + 2 i l ) ,,yv + ( x + 2 R) rny 

constants A and p (we use the notation 2 to distinguish it from 
the rake angle 2); y 3 is the horizontal displacement compo- 
nent of the eigenfunction; I\ and / 2 are energy integrals defined 
in Kanamori and Cipar [1974]. These functions need to be 
evaluated at the source depth d s to compute the change in grav- 
itational potential at the surface. The traction component of the 
eigenfunction, y 2 , must vanish at the surface and remains small 
at shallow depths. We can thus approximate n y 2l = A n X l 
+2 p n y l , rts 0, and obtain n y l ( ss — A t . In turn, the follow- 
ing approximation will hold for shallow sources: 




(21 + 1 ) 

4tc„ct ; 2 [/i + /(/ + 1 [A] 
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[ 93 ] Alternatively and introducing Poisson’s ratio 

-y == ^ 

2^/1 -f- /I s ) 5 


n K 0 u: 


(2Z+1) 


4n n crf[I\ + /(/ + l)/ 2 ] ll -2v 


„Xi(d Sj 


(A3) 


Note in particular that under this approximation, n Koj be- 
comes proportional to the excitation coefficient n Noj for an 
explosive source [Okal, 1978]: 


V n 


32- 


2/i 


2 p. 1 +v 

N a j — . _ n N 0 J . 


1 -2v 


(A4) 


[ 94 ] In the PREM model, (1 + v)/(l — 2v) varies only from 
2.6 to 2.9 within the upper and lower crust and the upper 
mantle. As a result, (A4) implies that dilatation in the source 
region will directly influence the isotropic gravity response 
function F 0 (d s ). 
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